1 Introducing the R Spatial Ecosystem (30 min)

1.1 The basics

1.1.1 rgdal, sp, rgeos

1.1.2 sf

  • Les fonctionnalités de sp, rgeos et rgdal dans un package unique.

  • Manipulation plus aisée, objets plus simples

  • Auteur principal et maintainer : Edzer Pebesma (auteur de sp)

  • Compatible avec les syntaxes pipe et les opérateurs du tidyverse.

Format des objets spatiaux sf

format sf

Import de données

library(sf)
mtq <- st_read("data/mtq/martinique.shp")
Reading layer `martinique' from data source `/home/tim/Documents/prz/flo/data/mtq/martinique.shp' using driver `ESRI Shapefile'
Simple feature collection with 34 features and 23 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
epsg (SRID):    32620
proj4string:    +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs

Affichage de données

plot(st_geometry(mtq))

Extraire les centroïdes

mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add=TRUE, cex=1.2, col="red", pch=20)

Construire une matrice de distances

mat <- st_distance(x=mtq_c,y=mtq_c)
mat[1:5,1:5]
Units: m
          [,1]     [,2]      [,3]      [,4]      [,5]
[1,]     0.000 35297.56  3091.501 12131.617 17136.310
[2,] 35297.557     0.00 38332.602 25518.913 18605.249
[3,]  3091.501 38332.60     0.000 15094.702 20226.198
[4,] 12131.617 25518.91 15094.702     0.000  7177.011
[5,] 17136.310 18605.25 20226.198  7177.011     0.000

Agréger des polygones

mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2, border = "red")

Construire une zone tampon

mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2)
plot(st_geometry(mtq_b), add=T, lwd=2, border = "red")

Réaliser une intersection

m <- rbind(c(700015,1624212), c(700015,1641586), c(719127,1641586), 
           c(719127,1624212), c(700015,1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border="red", lwd=2, add=T)

mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col="red", border="green", add=T)

Construire des polygones de Voronoi
google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r) 

mtq_v <- st_voronoi(x = st_union(mtq_c))
mtq_v <- st_intersection(st_cast(mtq_v), st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c, join=st_intersects)
mtq_v <- st_cast(mtq_v, "MULTIPOLYGON")
plot(st_geometry(mtq_v), col='lightblue')

1.2 Mapping

  • static:
    • ggplot2 + ggspatial
    • tmap
    • cartography
  • dynamic/interactive
    • leaflet
    • mapview

1.3 Spatial modelling and MISC

CRAN Task View: Analysis of Spatial Data

2 Exploration of sirene & osm db for restaurants

See file data_prep.R for data extraction.

# spatial data management
library(sf)
# cartography
library(cartography)
# smooth density computation
library(spatstat)
library(raster)
library(maptools)
# interactive mapping
library(mapview)

2.1 Display points

# load data
adm <- readRDS("data/iris_p13.rds")
feat_sir <- readRDS("data/sir_p13.rds")
feat_osm <- readRDS("data/osm_p13.rds")

nrow(feat_sir)
[1] 1066
nrow(feat_osm)
[1] 585
par(mar = c(0,0,1.2,0), mfrow = c(1,2))

plot(st_geometry(adm), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
plot(st_geometry(feat_sir), col = "red", pch = 20, add=T, cex = 0.5)
layoutLayer(title = "SIR", sources="", author="", scale = NULL, frame=FALSE)
legend(x = "topright", legend = "= 1 restaurant", pch = 20, pt.cex = 0.5, cex = 0.7, col = "red")

plot(st_geometry(adm), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
plot(st_geometry(feat_osm), col = "red", pch = 20, add=T, cex = 0.5)
layoutLayer(title = "OSM", sources="", author="", scale = .5, north = T, frame=FALSE)

2.2 Count points in admin units

inter_osm <- st_intersects(adm, feat_osm)
inter_sir <- st_intersects(adm, feat_sir)

adm <- st_sf(adm[,1:2,drop = T], 
             n_osm = sapply(X = inter_osm, FUN = length), 
             n_sir = sapply(X = inter_sir, FUN = length), 
             geometry = st_geometry(adm))

par(mar = c(0,0,1.2,0), mfrow = c(1,2))
plot(st_geometry(adm), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_sir", inches = 0.1)
layoutLayer(title = "SIR", sources="", author="", scale = NULL, frame=FALSE)

plot(st_geometry(adm), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_osm", inches = 0.1)
layoutLayer(title = "OSM", sources="", author="", scale = .5, north = T, frame=FALSE)

# plot(adm$n_sir, adm$n_osm, asp = T, xlim = c(0,55), pch = 21, bg = "red")
# text(adm$n_sir, adm$n_osm, labels = adm$CODE_IRIS, cex = 0.5, pos = 2)
# x <-lm(formula = adm$n_osm~adm$n_sir)
# abline(a = x$coefficients )

2.3 Interactive visualisation

Explore datasets interactively

The default is not bad:

mapview(feat_sir) + mapview(feat_osm, col.regions = "red")

mapview is hightly customizable:

mapview(feat_sir, map.types = "OpenStreetMap", col.regions = "#940000", 
        label = paste(feat_sir$L2_NORMALISEE, feat_sir$NOMEN_LONG, sep=" - "),
        color = "white", legend = TRUE, layer.name = "SIRENE", 
        homebutton = FALSE, lwd = 0.5) +
  mapview(feat_osm, col.regions = "#000094", color = "white", legend = TRUE, 
          label = feat_osm$name,
          lwd = 0.5, layer.name = "OSM",  homebutton = FALSE)

2.4 count points in regular grid

fufun2 <- function(feat,adm, cs = 500){
  grid <- st_make_grid(x = adm, cellsize = cs, what = "polygons")
  x <- st_intersects(grid, feat)
  gg <- st_sf(n = sapply(X = x, FUN = length), grid)
  plot(adm$geometry, col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
  choroLayer(gg[gg$n>0,], var = "n", border= NA, breaks = c(1,4,6,9,12,14,17,19,22), 
             col = carto.pal("wine.pal", 8), add=TRUE)
}
par(mar = c(0,0,1.2,0), mfrow = c(2,3))
fufun2(feat_sir, adm, 50)
fufun2(feat_sir, adm, 100)
fufun2(feat_sir, adm, 150)
fufun2(feat_osm, adm, 50)
fufun2(feat_osm, adm, 100)
fufun2(feat_osm, adm, 150)

2.5 smooth density analysis

Define a function that compute smoothed restaurant density and map it.

fufun <- function(feat, adm, title, sigma = 100,res = 50){
  bb <- as(feat, "Spatial")
  bbowin <- as.owin(as(adm, "Spatial"))
  pts <- coordinates(bb)
  p <- ppp(pts[,1], pts[,2], window=bbowin)
  ds <- density.ppp(p, sigma = sigma, eps = res)
  rasdens <- raster(ds) * 1000 * 1000
  proj4string(rasdens) <- "+init=epsg:2154"
  rasdens <- rasdens+1
  bks <- getBreaks(values(rasdens), nclass = 12, method = "equal")
  cols <- mapview::mapviewGetOption("raster.palette")(length(bks)-1)
  plot(adm$geometry, col = NA, border = NA, main="", bg = "#FBEDDA")
  plot(rasdens, breaks = bks, col=cols, add = T,legend=F)
  legendChoro(pos = "topright",cex = 0.7, title.cex = 0.7,
              title.txt = "resto per km2)",
              breaks = bks-1, nodata = FALSE,values.rnd = 0,
              col = cols)
  layoutLayer(title = paste0(title), scale = .5,
              tabtitle = TRUE, frame = FALSE,
              author = "Map data © OpenStreetMap contributors, under CC BY SA. SIRENE 2018", 
              sources = "cartography 2.0.2") 
  return(invisible(rasdens))
}
par(mar = c(0,0,1.2,0), mfrow = c(2,3))
fufun(feat_sir, adm, "SIR 50", 50, 20)
fufun(feat_sir, adm, "SIR 100", 100, 20)
fufun(feat_sir, adm, "SIR 250", 250, 20)
fufun(feat_osm, adm, "OSM 50", 50, 20)
fufun(feat_osm, adm, "OSM 100", 100, 20)
fufun(feat_osm, adm, "OSM 250", 250, 20)

reproducibility

sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 9 (stretch)

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.7.0
LAPACK: /usr/lib/lapack/liblapack.so.3.7.0

locale:
 [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C               LC_TIME=fr_FR.UTF-8       
 [4] LC_COLLATE=fr_FR.UTF-8     LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] mapview_2.4.0       maptools_0.9-2      raster_2.6-7        sp_1.3-1           
 [5] spatstat_1.56-0     rpart_4.1-13        nlme_3.1-137        spatstat.data_1.3-1
 [9] cartography_2.1.1   sf_0.6-3            knitr_1.20         

loaded via a namespace (and not attached):
 [1] lattice_0.20-35      colorspace_1.3-2     spatstat.utils_1.8-2 viridisLite_0.3.0   
 [5] stats4_3.5.1         htmltools_0.3.6      yaml_2.1.19          mgcv_1.8-24         
 [9] base64enc_0.1-3      R.oo_1.22.0          e1071_1.6-8          later_0.7.3         
[13] foreign_0.8-70       DBI_1.0.0            R.utils_2.6.0        RColorBrewer_1.1-2  
[17] plyr_1.8.4           foreach_1.4.4        stringr_1.3.1        rgeos_0.3-28        
[21] munsell_0.5.0        R.methodsS3_1.7.1    htmlwidgets_1.2      codetools_0.2-15    
[25] evaluate_0.10.1      httpuv_1.4.4.2       crosstalk_1.0.0      class_7.3-14        
[29] gdalUtils_2.0.1.14   Rcpp_0.12.17         xtable_1.8-2         tensor_1.5          
[33] scales_0.5.0         satellite_1.0.1      backports_1.1.2      promises_1.0.1      
[37] classInt_0.2-3       jsonlite_1.5         webshot_0.5.0        leaflet_2.0.1       
[41] abind_1.4-5          mime_0.5             deldir_0.1-15        png_0.1-7           
[45] digest_0.6.15        stringi_1.2.3        shiny_1.1.0          polyclip_1.9-0      
[49] grid_3.5.1           rprojroot_1.3-2      rgdal_1.3-3          tools_3.5.1         
[53] magrittr_1.5         goftest_1.1-1        Matrix_1.2-14        spData_0.2.9.0      
[57] rmarkdown_1.10       iterators_1.0.10     R6_2.2.2             units_0.6-0         
[61] compiler_3.5.1